library(lattice)
library(maps)

## If using fake Medicare data
# load(file="/Users/coryzigler/Dropbox/Research/HEI Accountability Grant/PM10 Accountability/Data/pm10_withps.Rda")

## If using real Medicare data
# load(file="/Users/coryzigler/Data/Medicare/Linked AQS-Medicare/HEI PM10 Accountability/pm10_withps.Rda")


#meanvars = c(paste("MEAN", 1987:1989, sep="."), paste("MEAN",1993:2011, sep="."))
whichyears = 1990:2001
whichvar = "Arithmetic.Mean" ##"X95th.Percentile" ##Arithmetic.Mean"
meanvars = c(paste(whichvar, whichyears, sep="."))
datwide = dat[, c("Monitor", "a", meanvars)]
names(datwide) = c("Monitor", "a", paste(whichvar, whichyears, sep="_"))
meanvars_r = c(paste(whichvar, whichyears, sep="_"))

datlong = reshape(datwide, direction="long", varying = meanvars_r, sep="_")
datlong = datlong[order(datlong$Monitor, datlong$time),]
datlong$MEAN = as.numeric(datlong[, whichvar])


a0means=with(subset(datlong, a==0), tapply(MEAN, time, mean, na.rm=TRUE))
a1means=with(subset(datlong, a==1), tapply(MEAN, time, mean, na.rm=TRUE))
plot(whichyears, a0means, ylim=range(c(a0means, a1means)), type="n", xlab="Year", ylab="PM10")
lines(whichyears, a0means, col=1, lwd=2)
lines(whichyears, a1means, col=2, lwd=2)


outdir = "/Users/coryzigler/Dropbox/Research/HEI Accountability Grant/Final Report/Report/Figures/"


#pdf(file = paste(outdir, "trackpoll1990_2001.pdf", sep=""))

lims = range(datlong$MEAN, na.rm=TRUE) 
linecolors = c(rgb(.5,.5,.5,.5), rgb(1,.5,.5,.5))
plot(range(whichyears), lims, type="n", xlab="", ylab="Annual Average Ambient PM10", axes=FALSE)
for (i in 1:dim(datwide)[[1]]){
 lines(whichyears, c(datwide[i, paste(whichvar, whichyears, sep="_")]), col=linecolors[datwide[i,"a"]+1], lwd=2)
}
lines(whichyears, a0means, lwd=2)
lines(whichyears, a1means, lwd=2, col=2)

axis(side=2, at=seq(5,105,5))
axis(side=1, at=whichyears, labels=whichyears)
mtext(side=1, line=2, at=1990, "Baseline")
mtext(side=1, line=2, at=2001, "Follow Up")

lines(whichyears, rep(50, length(whichyears)), lwd=2, lty=2, col=gray(.5))
#lines(rep(2001, 2), lims, lty=3, lwd=2, col=gray(.5))

legend("topleft", lty=c(1,1,2,3), lwd=c(2,2,2,2), col=c(1,2,gray(.5),gray(.5)), c("Attainment", "Nonattainment", "1987 NAAQS"), bty="n")

##add number of nonattainment monitors to plot
# Regulatory Data
attaindat=read.csv('/Users/coryzigler/Dropbox/Research/HEI Accountability Grant/Data/EPA Nonattainment Table.csv')
attaindat$FIPS=paste(formatC(attaindat$fips_state, width=2, flag="0"), formatC(attaindat$fips_cnty, width=3, flag="0"), sep="")
attaindat=subset(attaindat, pollutant=='PM-10')

countregs = merge(dat, attaindat, by="FIPS", all.x=TRUE)
regvars = paste("pw_", whichyears, sep="")
nreg = rep(NA, length(regvars))
for (i in 1:length(regvars))
	nreg[i] = sum(countregs[, regvars[i]] %in% c("P", "W"))

par(new=TRUE)
plot(whichyears, nreg, col=1, pch=16, axes=FALSE, xlab="", ylab="", type="n")
text(whichyears, nreg, nreg, cex=.65)
#axis(side=4, at = seq(220, 270, 5), labels= seq(220, 270, 5))

#dev.off()